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When  a  fine  powder  is  dumped  into  a  silo,  the  gas  trapped  by  the  particles  will  slowly  escape 
by  diffusing  through  the  material.  The  corresponding  uneven  gas  pressure  distribution 
creates  a  body  force  that  is  taken  into  account  through  Darcy's  law.  By  using  spatial 
averaging,  the  formulation,  even  though  essentially  one-dimensional  in  space,  includes 
effects  due  the  geometry  of  the  container.  An  efficient  and  robust  numerical  scheme  based 
on  a  DAE  formulation  is  proposed  and  implemented.  Various  computational  results  are 
presented  and  discussed  to  establish  the  validity  of  the  approach. 


I.  INTRODUCTION 

This  paper  deals  with  various  problems  related  to  the  simulation  of  aerated  powder 
consolidation.  This  kind  of  phenomena  is  routinely  encountered  during  the  handling  of 
fine  powders  in  countless  applications.  Typically  a  powder  is  stored  in  a  bunker  or  silo,  see 
Figure  1.  During  filling,  air  gets  trapped  in  the  material  leading,  in  some  cases,  to  partial 
fluidization  [8],  [14]  and  noticeable  changes  of  the  mechanical  properties.  Over  time  the 
excess  air  diffuses  through  the  powder  and  eventually  escapes  through  the  top  surface.  This 
paper  aims  to  find  the  length  of  time  a  given  material  takes  to  consolidate.  With  respect 
to  applications,  the  ultimate  goal  is  to  be  able  to  predict  and  consequently  avoid  flooding, 
i.e.,  the  sudden  discharge  from  a  hopper  of  a  fine  powder  at  a  much  greater  rate  than  that 
of  the  flow  of  ordinary  granular  materials.  Further  comments  on  the  connection  between 
flooding  and  deaeration  can  be  found  in  [12]. 

Early  work  by  Janssen  [6]  analyzed  the  behavior  of  a  column  of  granular  material  in  a 
container.  Forces  due  to  the  gradient  of  gas  pressure  were  neglected,  and  the  analysis  was 
restricted  to  vertical  cylinders.  Some  other  models  are  based  on  drastic  simplifications, 
such  as  considering  a  constant  vertical  stress  [9]. 

The  present  model  is  derived  from  basic  conservation  principles  of  mass  and  momentum. 
The  forces  resulting  from  nonuniform  gas  pressure  are  taken  into  account  through  Darcy’s 
law.  The  modeling  assumptions  roughly  follow  that  of  [8],  see  also  [4],  [12].  In  those 
publications,  the  effects  of  the  geometry  of  the  container  are  neglected,  in  effect,  treating 
only  the  case  of  cylindrical  bunkers.  In  [4],  that  case  is  mathematically  analyzed,  and  a 
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FIG.  1.  Geometry  and  coordinate  systems  for  the  vertical  conical  bunker.  The  height  of  the  column  of 
powder  of  time  t  is  denoted  H(t). 


robust  numerical  method  is  proposed  and  implemented.  One  of  the  contributions  of  this 
paper  is  the  extension  of  those  models  to  general  axisymmetric  domains.  Although  the 
modification  may  appear  as  a  technical  detail,  it  has  profound  ramifications  with  respect  to 
the  structure  of  the  system,  see  remarks  at  the  end  of  §2. 

Apart  from  using  specific  physical  constitutive  laws,  the  main  restrictive  assumption 
consists  of  neglecting  the  fluctuations  in  the  horizontal  directions,  allowing  for  an  essentially 
one  dimensional  in  space  formulation.  In  other  words,  the  problem  is  described  exclusively 
in  terms  of  quantities  that  have  been  averaged  in  the  horizontal  direction.  It  should  be 
noticed  that  in  many  situations  the  use  of  quasi  one-dimensional  consolidation  models  can 
only  be  viewed  as  a  first  step,  and  full  multidimensional  approaches  should  be  considered 
instead.  We  refer  to  [  1 6]  and  [10],  Chap.5,  for  remarks  about  the  limitations  of  such  models. 
More  importantly  however,  even  a  simplified  model  such  as  the  one  considered  here  is  very 
sensitive  to  the  values  of  various  material  coefficients  such  as  compressibility.  This  is 
clearly  illustrated  by  our  numerical  results,  see  e.g.  Figure  5.  Those  material  coefficients 
are  typically  hard  to  measure  in  an  accurate  way. 

The  paper  is  organized  as  follows.  The  model  is  derived  in  Section  2.  The  resulting 
system  is  nonlinear  and  strongly  coupled.  It  consists  of  an  essentially  parabolic  PDE,  an 
ODE  and  an  integral  equation.  The  three  main  unknowns  are  the  average  vertical  stress,  the 
average  gas  pressure  and  the  height  of  the  powder  in  the  container.  Through  an  appropriate 
transformation,  the  calculations  are  performed  in  a  fixed  reference  computational  domain. 
A  discretization  is  proposed  in  Section  3.  The  spatial  discretization  is  second  order  accurate 
and  uses  a  combination  of  centered  Finite  Differences  and  BDF  [  1  ],  [5].  The  semi  discretized 
in  space  system  corresponds  to  a  semi-explicit  index  2  Differential  Algebraic  Equation 
(DAE).  The  time  discretization  is  done  by  a  linearly  implicit  Euler  discretization,  which, 
although  only  first  order,  is  the  simplest  acceptable  numerical  approach  for  the  above  type 
of  DAEs.  Our  numerical  experiments  show  it  to  be  both  robust  and  efficient  in  the  present 
context.  Computational  results  are  presented  and  discussed  Section  4.  Finally,  Section  5  is 
devoted  to  concluding  remarks. 


2.  THE  MODEL 

We  consider  general  axisymmetric  hoppers  as  illustrated  in  Figure  1 .  In  order  to  simplify 
the  problem,  a  pseudo  one-dimensional  formulation  is  derived  by  averaging  all  quantities 
on  horizontal  cross-sections.  The  height  of  the  powder  in  the  container  at  a  generic  time 
t  can  thus  be  described  by  a  function  H(t).  Using  cylindrical  coordinates,  the  domain 


POWDER  CONSOLIDATION 


3 


occupied  by  the  powder  at  time  t  is  given  by 

{(r,  0,  z);  0  <  z  <  H(t),  0<6><27r,0<r<  R(z )}, 

where  R(z)  stands  for  the  radius  of  the  hopper  at  height  z. 

For  an  arbitrary  function  /,  we  define 

i  rR(z)  r27r 

L  f^rdrdd’ 

where  f(z)  is  the  averaged  value  of  the  function  /  at  height  z  in  the  hopper.  By  axisymmetry, 
the  functions  considered  here  do  not  depend  on  the  angular  variable  6 ,  and  thus  with  a  slight 
abuse  of  notation 


Invoking  axisymmetry,  the  stress  tensor  has  the  form 

(Xrr  0  uTZ  1 

T  =  0  aee  0 

_  ®rz  0  ® zz  _ 

Consider  an  infinitesimal  slice  of  material  of  height  Sz.  The  forces  acting  on  such  a  slice 
are  averaged  and  summed.  The  average  vertical  stress  is  clearly  given  by  azz.  We  denote 
by  7  and  p  the  bulk  density  and  gas  pressure  respectively,  and  by  7  and  p,  their  respective 
average  values.  The  various  forces  are 

•  weight  of  solid:  -77r  R2(z)Sz; 

•  if  fw  is  the  average  wall  shear  stress  and  aw  the  average  normal  stress  on  the  wall, 
there  are  upward  forces  of  2  7r  R{z)  Sz  fw  and  2  7 r  R(z)  R'(z)  8zaw ; 

•  upward  pressure  due  to  the  wall:  2  n  R(z)  R'(z)  Szp(z ,  t); 

•  pressure  at  bottom:  p(z,t)  7r  R2(z),  and  top:  —  (p(z,  t)+8p)  7r  R2(z+Sz);  this  creates 
a  force  — 7r  R(z)2  Sp  —  2  tt  R(z)  R'(z)  p(z ,  t)  Sz; 

•  average  vertical  stress  at  bottom:  ozz(z)^  and  top:  —{<Jzz(z)  +  8crzz);  (compressive 
stresses  are  taken  as  positive  for  granular  material);  this  creates  a  force  — 7r  R2(z)  Sgzz  — 
2ix  R(z)  Rf(z)  azz  Sz. 

The  resulting  balance  of  forces  equation  is 

0  o 

dzGZz  +  dzP  H — R(z)~  ^ zz  ~  R(z )  +  7  ”  0. 

On  the  wall,  the  law  of  sliding  friction  applies 

Tyj  =  pW&Wl 

where  fiw  is  the  coefficient  of  wall  friction.  The  average  stress  tensor  T  is  diagonal; 
in  other  words,  both  orr  and  dzz  are  principal  while  arz  —  0.  Notice  that  aw(z)  = 
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cos2  0(z)arr(z)  4-  sin2  6(z)gzz(z)  where  tan  0(z)  —  R'(z).  The  ratio  of  average  vertical 
stress  ozz  to  average  wall  stress  gw 


—  1  V-/ 

&ZZ 

is  taken  as  depending  on  the  geometry  only,  see  §4  for  more  details.  The  above  assumption 
is  the  pendant  to  Janssen’s  analysis,  [6],  [8],  [10],  which  is  routinely  used  in  vertical 
bunkers.  In  (3),  k  depends  on  the  geometry  of  the  container  and  has  to  be  determined 
through  experiments  [8],  [13]  and  §4.  Equations  (2)  and  (3)  then  yield 


dzG  zz  +  dzp  4- 


2R'(z)  _  2 

R(z)  °zz  R{z) 


(pw  4-  R'(z))  kazz-\-  7  =  0. 


(4) 


Let  T  be  the  density  of  the  solid  particles,  which  is  assumed  constant.  The  gas  density, 
denoted  by  p,  is  an  unknown  function  of  time  and  position.  These  two  quantities  are  linked 
through  the  bulk  density  7 


7  =  fsT  4-  (1  —  fs)p,  (5) 

where  fs  is  the  volume  fraction  occupied  by  the  solid.  Generally  F  is  at  least  three  orders 
of  magnitude  larger  that  p,  and  thus 


/.  «  f-  (6) 

The  average  bulk  density  is  considered  as  a  function  of  the  average  major  consolidating 
stress,  here  aZZi  i.e.,  7  =  7( dzz ).  Various  such  relations  have  been  proposed,  see  [10], 
§6.2,  for  a  review.  Those  models  typically  make  sense  for  only  a  limited  range  of  values  of 
gzz ;.  Following  [7],  we  assume 


7  =  7m(l  +  — 4,  (7) 

Gm 

where  0  <  (3  <  1  is  the  coefficient  of  compressibility  of  the  material,  >  0  and  Gm  >  0 
are  material  constants.  We  refer  to  [3]  and  [2]  for  respectively  theoretical  and  experimental 
investigations  of  the  precise  nature  of  the  bulk  density/stress  relation. 

Assuming  powder  is  neither  entering  nor  leaving  the  system,  the  total  mass  M  of  the 
solid  is  conserved,  leading  to 

pH (*)  M 

/  R2(z)  j(G(z,t))  dz  =  —  =  M,  t  >  0.  (8) 

Jo  K 

Applying  the  continuity  equation  to  both  gas  and  solid  phases,  we  get 

dtl  4-  V  •  (7  us)  =  0 

dt  ((1  -  £)  p)  +  v  •  ((1  -  2)  p  Ug)  =  0, 


where  us  and  ug  are  the  velocities  of  the  solid  and  the  gas,  respectively. 


(9) 

(10) 
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Using  axial  symmetry  and  applying  the  averaging  operator  (1)  to  (9)  yields 

2  2  fR(z)  /  \ 

dt 7  +  z )  UsAR(z)> z )  +  Jo  dz  z)  u*Ari  z)j  rdr  =  0. 

Further,  assuming  the  grains  in  contact  with  the  wall  to  move  tangentially  with  respect  to 
it,  we  observe  us^r(R(z),  z)  =  R!(z)  us^z(R(z),  z).  Elementary  manipulations  then  leads 
to 

df 7  +  dz  ( R 2 (z) lulA)  =  0. 

Neglecting  fluctuations  in  the  radial  direction  gives  7  us,z  «  7  u^.  Similar  principles  can 
be  applied  to  local  conservation  of  gas.  Those  conservation  laws  then  read 


dt7  +  W{z)dz  (jR2^7Ws’z1  =  °’ 


(ii) 


dt(K{l-l)p)  +  1^)dz(R\z)p{l-l)ug,z)  =  0.  (12) 

In  addition  to  (3)  and  (7),  two  additional  constitutive  equations  are  considered.  First,  the 
gas  is  assumed  to  be  ideal  and  isothermal 


P  =  Po 
P  Po ’ 


(13) 


where  po  and  po  are  constant  reference  values.  Second,  pressure  gradient  and  velocities 
are  related  through  Darcy’s  law,  which  reads  here 

&g,z  “  Us,z  =  “^(7)  SzP,  (14) 

where  K  is  the  permeability,  taken  as  a  function  of  the  average  bulk  density.  We  take  [8] 


*«)=*»  (7"' 


(15) 


where  Kq  and  70  are  reference  values  and  a  is  a  positive  constant.  The  parameters  (3,  <xm, 
7m,  a,  K0  and  70  appearing  in  (7)  and  (15)  have  to  be  determined  experimentally. 

We  now  eliminate  velocities  from  the  system.  Putting  together  (14)  and  (1 1)  yields,  after 
integration  between  0  and  a  generic  point  z 

R2(z)da(z)dz  +  R2(z)Az)(u9Az)  +  K(l(z))dzp(z)j 

-  -R2 (0)7(0) (%2(0)  +  tf(7(0))ey>(0))  =  0. 

Assuming  no  gas  leaves  the  (closed)  container  at  2  =  0,  natural  boundary  conditions  at 
the  bottom  are  then  u9iZ( 0,  t)  —  0  and  dzp( 0,  t)  =  0.  The  previous  equation  can  be  solved 
for  u9tZ(z),  leading  to 

Ug,z(Z)  =  no2?"  T  [  R2W  dt^  dz  ~  KA(Z))  9zp(z). 

y  7 (z)R2(z)  Jo 
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Using  (13)  and  the  above  expression  for  ugjZ(z)9  (12)  can  be  rewritten  in  terms  of  pressure 
p  eliminating  p  from  the  problem. 

The  full  system  of  equations  can  now  be  considered.  The  three  unknowns  of  the 
problem  are  the  averaged  gas  pressure  p ,  the  averaged  vertical  stress  azz  and  the  height 
of  the  powder  H(t).  We  obtain  the  following  system  of  integrodifferential  equations  in 
{(z,  £);0  <  t,  0  <  z  <  H(t)} 

o  -  m 

-(1  -  1)  dz(pKdzp)  +  dzp, 

0  =  dzazz  +  dzp+2^^azz-~-^(iJ,w  +  R'(z))kazz  +  j,  (17) 
tW) 

M  =  /  R2(z)i(azz(z,t))  dz.  (18) 

Jo 

The  system  is  completed  by  initial  and  boundary  conditions 

p(-iO)  =  Po, 

azz(H(t),t)  =  0,  dzpiOjt)  =  0,  =patm, 


where  patm  is  the  value  of  atmospheric  pressure.  The  boundary  condition  at  z  =  0 
corresponds  to  a  solid  bottomed  container.  This  can  be  easily  modified  to  include  the  case 
when  some  material  exits  through  an  outlet  at  the  bottom.  A  mathematical  analysis  of  the 
above  problem  in  the  case  of  a  cylindrical  container,  i.e.,  R'(z)  =  0,  can  be  found  in  [4]. 

A  useful  tool  for  the  present  kind  of  equations  consists  in  mapping  the  entire  problem 
into  a  fixed  spatial  domain.  Upon  discretization,  this  type  of  formulation  leads  to  better 
stability  properties  than  formulations  where  the  grid  moves  with  the  material,  such  as  [8]. 
Let 

y  =  z/H(t)  p(y,  t)  =  p(z,  t)  (r(y,t)  =  azz(z,t). 

Based  on  these  relations,  the  previous  system  of  equations  (16),  (17)  and  (18)  yields 


0  -  (1- f)  (ftp -9 (va*-»vvf$ 

-  Tmrnd>  K  “  ?>)  f R2{Him  (ya<’ "  S7'fl'£’§i) 

“  Tlkn1'  ~  f  +  rnw,'/i,’c'd’p' 

0  =  dy(T  +  dyP  +  2 


(19) 


xymbm.  _  + *w)  + mt)l,  (20) 


- a  — 


R(H(t)y)  ”  R(H(t)y) 

M  =  H(t)  f  R2(H(t)y)j(a(y,t))dy, 

Jo 

where  we  have  dropped  the  bar  notation.  Initial  and  boundary  conditions  become 


(21) 


p(y,  0)  =  po(y),  y  e  (0,  l), 

CT(i,t)  =  0,  dyp{0,t)  =  0,  p(i,t)=patm, 


(22) 

(23) 
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The  change  from  R(z)  =  cst  to  a  general  function  R(z)  is  far  from  being  benign.  If  #(2)  = 
iTo  =  cst,  (21)  can  be  solved  for  H(t):  H(t)  =  (f*  7 (cr(y,t))  dy^J  .  Taking  the 

derivative,  one  can  also  obtain  an  explicit  expression  for  H'(t).  Those  expressions  can  be 
then  substituted  in  (19)  and  (20)  to  fully  eliminate  H  from  the  system  [4].  For  general 
functions  R  however,  (21)  is  a  nonlinear  equation  for  H(t),  which  cannot  be  expected  to 
always  yield  explicit  solutions.  To  maintain  generality,  H(t)  is  not  eliminated  from  the 
system  before  discretization  but  is  rather  computed  as  an  unknown  along  with  stress  and 
pressure. 


3.  DISCRETIZATION 

These  equations  are  discretized  and  solved  in  the  fixed  rectangular  domain  (0, 1)  x  (0,  T). 
Set  Ay  =  1/N  and  let  ^  =  (i  —  1)A y  i  =  1, . . . ,  N  +  1.  The  semidiscretized  variables 
are  pressure  P(t)  =  [Pi(t), . . . ,  Pjv(t) ]  and  stress  £(£)  =  [£i(£), . . . ,  EjvOfc)].  Note  that 
the  remaining  unknown,  iJ,  depends  only  on  time.  The  boundary  conditions  at  y  =  1  read 
PN+i  =  Patm  and  Ejv+i  =0.  A  second  order  spatial  discretization  for  both  variables  is 
used.  The  discretized  operators  are 

Cp :  N  x  N  matrix  corresponding  to  a  second  order  centered  discretization  of  dy  with 
pressure  boundary  conditions  (at  y  —  0  and  y  —  1), 

D:  N  x  N  matrix  corresponding  to  a  second  order  centered  discretization  of  dyy  with 
pressure  boundary  conditions  (at  y  =  0  and  y  =  1), 

Ca:  N  x  N  matrix  corresponding  to  a  second  order  backward  differentiation  formula 
(BDF  [1],  [5])  for  dy  with  stress  boundary  condition  (at  y  —  1). 

Although  more  sophisticated  discretization  methods,  such  as  BDF  based  algorithms  [1], 
should  be  considered  in  space  as  well  in  general,  the  diffusion  like  dynamics  of  the  present 
problem  are  sufficiently  simple  to  mitigate  this  point  of  view.  The  construction  of  Cp  and 
D  is  elementary.  For  Ca  we  take 


3  4-1  0  .  0 

0-3  4-1  0  ...  0 


0  ...  0-3  4-1  0  . 

0  .  0-3  4-1 

0  .  0-3  4 

0  .  0  -2_ 

In  addition,  the  following  discretized  integral  operators  are  introduced 

Aa(W)  =  UjWj  md(BA(W))i  =  ^2ujjWj, 

M  j= 1  j'=i 

where  u>  denotes  the  weights  of  the  quadrature.  Three  additional  variables  corresponding  to 
the  derivatives  of  the  three  main  three  unknowns  P, ).'  and  H  are  also  introduced,  U  =  dtP, 
V  =  d(E  and  L  =  H’.  The  spatially  discretized  problem  is 


dtP  =U 


(24) 
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o  =  (i  - 1)  (u  -  §y  cpp)  -  p  (y  v  -  y |  y  c„e)  (25) 

~w  ((^  "  f)CpP  ~  ^E)  -  l^^yyc.E)) 

1  *y  PK 

~(1  -  ^(P,  E)P  +  f^j'C^CpP  +  tel  =  /(P,  U,  V,  E,  H, L) 

O  TT  T>f  9  L7 

0  =CCTE  +  CpP  +  -^7yE  +  6c2--^Fy(^  +  P')  +  ^7  (26) 

=  G(P,U,V,Z,H,L ) 

dtE  =  V  (27) 

dtP  =  L  (28) 

0  =  P.4A  (P2(PY)7)  s  ft(P,  U,  V. ;  E,  H,  L),  (29) 


where  7(E)  and  7'(E)  are  vectors  and  Y  =  [j/i, . .  - ,  2/jv]-  The  two  vectors  hcl  and  bc2 
result  from  the  presence  of  boundary  conditions.  Also,  vector  to  vector  multiplication  in 
the  above  expressions  is  taken  to  be  componentwise.  The  problem  has  the  form 

dtP  =  U, 

0  =  /(P,  U,  V,  E,  H,  L), 

0  =  g(P,0, 0,  E,iT,0), 
dt  E  =  V 
dtH  =  L 

0  =  ft(0, 0, 0,  E,  H,  0) 


This  corresponds  to  an  semi-explicit  index  2  DAE  or  equivalently  a  fully  implicit  index  1 
DAE  [1].  A  linearly  Implicit  Euler  time  discretization  [5],  [1 1]  of  (24)-(29)  leads 


I  -A  tl  0  0  0  0  ‘ 

pn+l  _  pn 

■  jjn  - 

-A t/r  —  At/2  -Af/3n  -A tft  -At  ft  -A t/6n 

un+ 1  -  un 

r 

—Atgi  -A tfi  -AtgS  -At  iff  -At  gg  -At  gg 

yn+1  _  yn 

=  At 

9n 

O 

O 

1 

> 

?*+• 

O 

O 

Sn+1  _  Sn 

yn 

0  0  0  0  1  -At 

fjn+1  _  Hn 

Ln 

-At  kg  -Ath%  -Athg  -A thg  -At  ft?  -A thg  _ 

Ln+ 1  -  Ln 

.  h"  . 

where  the  subscripts  denote  derivatives  with  respect  to  the  corresponding  variables.  The 
first  partial  derivatives  of  /  and  g  are  evaluated  through  finite  differences. 


4.  COMPUTATIONAL  RESULTS 

Two  cases  with  different  geometries  are  considered.  The  first  corresponds  to  a  cylinder 
and  the  second  to  a  bunker  consisting  of  a  cylinder  for  the  upper  portion  and  a  cone  for  the 
lower  portion.  The  data  is  initialized  as  follows. 

The  values  of  various  parameters  in  the  model  have  been  taken  as  to  represent  a  realistic 
situation 


An  =  -25 
7m  =  60  lbs/ft3 
am  =  13  lbs/ft2 
7o  =  80  lbs/ft3 


T  =  200  lbs/ft3 

K0  =  10"4  f^lbs-V1 

G  —  4 

=  tan(7r/9) 
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The  atmospheric  pressure  is  patm  =  2116.8  lbs/ft2 . 

The  determination  of  the  ratio  k  between  wall  stress  and  vertical  stress,  see  (3),  is  more 
delicate.  In  the  case  of  a  perfectly  plastic  flow,  such  a  coefficient  could  be  derived  from 
a  plasticity  model,  Mohr-Coulomb  for  instance  [10].  Although  this  is  often  done  as  part 
of  an  analysis  a  la  Janssen,  such  an  approach  is  not  consistent  with  the  mechanical  states 
considered  here  since  one  cannot  assume  the  material  to  have  reached  yield.  The  value  of 
the  ratio  k  does,  however,  depend  on  the  geometry.  For  plastic  granular  flows,  the  stresses 
are  usually  believed  to  be  close  to  an  active  state  in  a  vertical  cylindrical  hopper,  while 
they  tend  to  be  in  a  passive  state  in  converging  conical  hoppers  for  instance  [10],  [15]. 
Experimental  evidence  [13]  seems  to  point  to  the  fact  that,  even  here,  a  similar  behavior 
can  be  observed.  More  precisely,  the  value  of  k  tends  to  be  much  lower  in  a  converging 
conical  section  than  in  a  vertical  cylindrical  one.  The  values  below  are  taken  to  reflect  this. 

Let  us  consider  the  cylinder  with  radius  R  =  0.5ft  (and  thus  R'  =  0)  with  an  initial 
height  of  powder  iif(O)  =  5ft.  The  ratio  k  is  taken  as  k  =  1/3.  An  initial  pressure  field, 
po,  consistent  with  the  boundary  conditions  in  (23)  is  given  as 

Po(z)  =  Patm  (l  +  100  (x  ”  (#(0))  ))  • 

In  the  calculations  below,  N  —  200  and  NT  =  T/At  —  500. 

For  the  second  case,  the  radius  of  the  cylinder  is  R  =  0.5  ft.  The  radius  at  the  bottom 
of  the  bunker  is  0.4ft,  and  the  height  of  the  conical  part  is  1ft.  Thus,  R'  =  0.1.  In  order 
to  have  approximately  the  same  mass  of  material  in  this  hopper  as  the  first  case,  the  total 
initial  height  of  the  powder  is  H( 0)  =  5.459ft.  The  considered  values  of  k  are  1/3  in  the 
cylindrical  part,  and  k  =  3  in  the  conical  part.  The  same  initial  pressure  field  is  taken.  In 
addition,  N  =  200  and  NT  —  T/At  =  500  as  before. 

Figure  2  gives  a  typical  illustration  of  the  behavior  of  the  problem  in  a  cylinder.  Note  that 
Figure  2  is  related  to  the  original  geometry  of  the  problem,  i.e.,  the  variable  is  z  not  y.  As 
expected,  a  settlement  of  the  powder  is  observed.  For  both  graphs,  the  top  line  corresponds 
to  the  position  of  the  top  of  the  powder  column,  i.e.,  the  free  boundary.  Further,  the  pressure 
field  is  found  to  asymptotically  converge  to  a  uniform  pressure  value  corresponding  to  patm, 
the  atmospheric  pressure,  i.e.,  equilibrium  of  the  pressure  is  established.  The  stress  field 
is  also  found  to  converge  to  a  stationary  distribution  which,  in  terms  of  the  original 
variable  2,  is  solution  to 


dzaa 


R(z) 


(lyj  kcToo  -j~  tC^oo) 

^oo(-^oo)  =  0, 


0  z  H0 o, 


where  Hoo  is  the  asymptotic  value  of  the  height  of  the  powder  column.  Figure  4,  left, 
illustrates  the  convergence  to  the  asymptotic  stationary  state  after  2000s. 

Figure  3  gives  a  typical  illustration  of  the  behavior  of  the  problem  in  a  bunker  consisting 
of  a  cylinder  as  the  upper  part  and  a  cone  in  the  lower  part.  Again  settlement  of  the  powder 
is  observed  by  considering  the  top  line  of  the  graphs.  As  in  the  previous  case,  the  pressure 
field  asymptotically  converges  to  a  uniform  pressure  value  corresponding  to  patm  •  The 
stress  field  also  converges  to  a  stationaiy  distribution  <7^  which,  in  terms  of  the  original 
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FIG.  2.  Calculated  pressure  and  stress  fields  in  the  cylindrical  geometry;  N  =  200,  NT  —  500,  T  —  1500. 
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FIG.  3.  Calculated  pressure  and  stress  fields  in  the  cylinder-on-a-cone  geometry;  N  =  200,  NT  =  500, 
T  =  1500s. 
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FIG.  4.  Calculated  stress  field  in  the  original  geometry;  left:  cylinder,  right: cylinder/cone  geometry. 


variable  z,  is  solution  to 


*  R' 

dZ(7 oo  +  2-— rCToo  -  • 


'R(z)  °°  R(z) 


(flw  ^R^kcToo  +  7(^00 )  =  o,  0<z<iJo, 
CToo(tfoo)  =  0? 


where  Hoo  is  the  asymptotic  value  of  the  height  of  the  powder  column.  Figure  4  illustrates 
the  convergence  to  the  asymptotic  stationary  state  after  2000s. 

To  verify  the  results  of  the  calculations,  an  experiment  must  be  designed  so  that  the 
parameters  for  the  model  can  be  determined.  Very  accurate  measurements  of  these  param¬ 
eters  are  needed  to  verify  the  calculated  results  with  experimental  ones.  Consider  Figure  5 
which  shows  the  half-life  of  the  over-pressure  at  the  bottom  of  a  cylindrical  bunker  versus 
the  compressibility  parameter  (3 .  The  half-lives  range  from  nearly  zero  to  approximately 
650  for  (3  ranging  from  zero  to  four.  It  can  be  shown  that  the  half-life  of  the  over-pressure 
is  also  sensitive  to  other  parameters  in  the  model  [12]. 


5.  CONCLUSIONS 

Existing  models  for  powder  consolidation  have  been  extended  to  general  axisymmetric 
domains  in  order  to  take  into  account  geometrical  effects.  In  this  process  and  as  with  those 
prior  models,  averaging  on  horizontal  cross  sections  plays  a  fundamental  simplifying  role. 
The  resulting  system  consists  of  an  essentially  parabolic  PDE,  and  an  ODE  in  space  and 
an  integral  equation.  In  spite  of  the  strong  nonlinear  coupling,  those  can  be  thought  of  as 
respectively  equations  for  the  gas  pressure,  the  vertical  stress  and  the  height  of  the  powder. 
A  fundamental  and  problematic  material  parameter  is  the  ratio  k  between  wall  stress  and 
vertical  stress. 

Unfortunately,  some  of  the  parameters  entering  the  model,  such  as  k  or  the  compressibil¬ 
ity  /?,  are  not  easy  to  measure  in  a  reliable  way.  To  make  matters  worse,  the  present  work 
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FIG.  5.  Calculated  half-life  of  over-pressure  versus  the  parameter  / 3 . 


illustrates  that  even  relatively  small  variations  in  those  coefficients  can  have  a  noticeable 
effect  on  the  solution,  see,  e.g.,  Figure  5.  This  may  partially  explain  the  relative  sparsity  of 
experimental  results  for  the  present  type  of  problems. 

The  main  contributions  of  the  paper  are  first  a  careful  generalization  of  earlier  models  in 
order  to  take  geometrical  effects  into  account  and  second  the  design  and  implementation 
of  a  robust  and  efficient  numerical  algorithm  for  the  corresponding  equations.  While  they 
extend  earlier  results  [4],  [8],  [12],  the  presented  computational  experiments  seem  to  be 
consistent  with  them. 

Various  aspects  of  the  present  model  are  questionable,  primarily  the  use  of  a  pseudo  one¬ 
dimensional  formulation  through  spatial  averaging.  Although  a  truly  multidimensional 
approach  appears  to  be  desirable,  its  precise  nature,  especially  in  term  of  constitutive 
assumptions,  is  unclear. 
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